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I.   INTRODUCTION 

Many  important  structural  problems  In  aerospace  engi- 
neering involve  components  that  are  in  contact  with  a  fluid. 
One  example  is  the  fuel  tank  located  in  a  wing  or  fuselage 
of  an  aircraft.   Another  example  is  the  fuel  tank  in  a 
liquid  fueled  rocket.   An  example  in  the  field  of  civil 
engineering  is  the  dam  and  reservoir  system. 

An  ongoing  research  program  at  the  Naval  Postgraduate 
School  is  concerned  with  aircraft  vulnerability  in  the  fuel 
tank  area.   This  program  is  currently  devoted  to  a  study  of 
the  Hydraulic  Ram  Problem.   "Hydraulic  Ram"  is  the  dynamic 
loading  of  fuel  tanks  when  impacted  by  bullets  or  warhead 
fragments.   This  loading  gives  the  appearance  of  a  simple 
pressure  pulse  with  a  resulting  deformation  or  rupture  of 
the  tank.   It  has  been  found,  however,  that  the  loading  and 
deformation  are  a  combination  of  a  number  of  events  called 
"hydraulic  ram  components."  [Ref.  1] 

As  part  of  the  theoretical  analysis  of  the  fuel  tank 
response  to  a  penetrating  projectile  two  computer  codes  were 
written  to  solve  for  the  response  of  the  fluid  and  the  walls 
of  a  idealized  two-dimensional  fuel  tank.   The  codes  will 
compute  the  pressure  throughout  the  fluid  and  will  calculate 
the  wall  response  to  this  pressure.   The  user  can  prescribe 
the  problem  to  be  analyzed,  includings but  not  limited  to: 


Fluid  surface  pressure  disturbances,  initial  wall  displace- 
ments, initial  wall  velocities,  pressure  distributions  through 
the  fluid,  etc. 

This  thesis  presents  the  method  selected  for  the  solution 
of  problems  involving  two-dimensional  fluid-structure  inter- 
action.  The  types  of  problems  which  can  be  solved  using  the 
codes  are  described.   Section  II  contains  the  theoretical 
description  of  the  problem,  and  the  following  two  sections 
outline  two  finite  difference  formulations  that  were  used 
to  obtain  solutions.   Section  V  contains  a  description  of 
each  code  and  describes  the  modifications  necessary  to 
change  the  boundary  and  initial  conditions.   Section  VI 
presents  the  solutions  obtained  with  the  codes,  including 
the  check  cases  and  several  example  programs,  and  discusses 
the  comparison  of  results  obtained  from  the  two  code  formula- 
tions.  Users'  instructions  and  program  listings  for  both 
formulations  are  included  in  Appendices  A  and  B. 


II.   DERIVATION  OF  GOVERNING  EQUATIONS 

A.   THE  WAVE  EQUATION  FOR  AN  INVISCID  FLUID 
1.   Equation  of  Motion 

Euler's  equation  of  motion  can  be  given  in  the  form 


P  §f  +  grad  p  =   0  (1) 


where  P  is  the  fluid  density,  p  is  the  fluid  hydrostatic 
pressure  (positive  in  compression),  q  is  the  velocity  vector 
and  t  is  time.   The  total  derivative  =jr-  is  given  by 


bit  =  !t  +  ufe+  vi?  +  wft  <2> 


where  u,  v,  w  are  the  components  of  the  velocity  in  the 

x,  y,  z  directions  respectively.   When  the  fluid  velocity  is 

sufficiently  small, 


2-  *  £-  (3) 

Dt   at  VJ; 


Superscript  dots  will  be  used  to  denote  partial  derivatives 
with  respect  to  time. 

2 .   Constitutive  Equation 

The  constitutive  equation  of  the  inviscid  fluid  is 
given  by 

p  =  -  K  div  q  (4) 

where  K  is  the  bulk  modulus. 


3.   Equation  of  Continuity 


ij^+  divq   =   0  (5) 


Applying  Eq .  (3)  we  obtain 


-  p  +  div  q   =   0  (6) 


Eliminating  div  q  between  Eqs .  (4)  and  (6)  gives 


K    p 


4 .   Potential  Equation 
Define  a  potential 


q  ^  grad  tj>  (8) 


Thus  Eq .  (1)  becomes 

p  grad  $  +  grad  p  =   0  (9) 

when  Eq .  (3)  is  applied.   Assuming  p  to  be  essentially 
constant  throughout  the  fluid,  Eq .  (9)  becomes 

grad  (p4>  +  p)   =   0  (10) 


or 


p$  +  p   =   constant  (11) 
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The  constant  term  is  set  to  zero  since  it  implies  a  base 
pressure  which  is  constant  throughout  the  domain. 

Differentiating  Eq .  (1)  with  respect  to  time  gives 

pq  +  grad  p   =   0  (12) 

Applying  Eqs .  (4)  and  (8) 

p  grad  <J)  +  grad  (-K  div  grad  <f>)   =   0      (13) 

Since  p  is  assumed  to  be  essentially  constant,  Eq .  (13) 
becomes 

grad(p<J)  -  K  div  grad  <j>)   =   0  (1*0 


or 


grad(p<{>  -  KV2<}>)      =      0  (15) 


Then 


2 
p<f>   -   KV   <{>     =      constant  (16) 


or 


p<j>   =   KV2<f>  (17) 


Recalling  that  the  speed  of  sound  in  a  fluid  can  be  deter- 
mined from  the  relationship 


c2  =  £  (18) 

P 
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we    obtain   the  wave   equation 


c2V2<{>     =      <J>  (19) 


B.  THE  GOVERNING  EQUATION  FOR  THE  STRUCTURE 
1.   Two-Dimensional  Idealization 

Consider  a  long  open  box  composed  of  five  uniform 
plates.   It  can  be  idealized  in  two  dimensions  as  a  structure 
composed  of  three  unit  width  beams  attached  at  right  angles 
as  shown  in  Fig.  1.   In  the  analysis  each  beam  is  considered 
separately,  with  its  own  set  of  independent  boundary  conditions 
In  all  of  the  examples  considered,  each  beam  was  assumed  to 
be  clamped  at  both  ends.' 

The  equation  for  the 'lateral  displacement  of  a  uniform 
beam  under  dynamic  load  is  [Ref.  2] 

EI  — jj-  w(x,t)  +  mw(x,t)  =  -  p  (20) 

9x 

where  E  is  Young's  modulus,  I  is  the  section  moment  of 
inertia,  m  is  the  mass  per  unit  length,  w(x,t)  is  the  lateral 
displacement  of  the  beam,  and  p  is  the  pressure  applied  to 
the  beam.   A  positive  pressure  is  in  the  opposite  direction 
as  w(x,t)  as  shown  in  Fig.  1  for  the  tank. 

C.  INTERFACE  CONDITIONS 

From  Eq .  (11)  the  pressure  on  the  beam  is  the  fluid 
pressure  given  by 
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P  =  -  P4>  (21) 

where  again  a  positive  pressure  is  in  the  opposite  direction 
as  w(x,t) . 

Another  interface  condition  is  the  requirement  that  the 
normal  velocity  of  the  fluid  potential  at  the  interface  is 
equal  to  the  velocity  of  the  wall,  i.e. 


w(x,t)  =  K   <J>  (22) 

8n 

where  — —  indicates  a  partial  derivative  normal  to  the  wall. 
8n 

Where  the  fluid  has  a  free  surface  and  the  pressure  is 

zero,  the  boundary  condition  is  given  by 


*  =  0  (23) 


When  an  excitation  pressure  is  applied  to  the  free  surface 


I   =  -  E.  (2*0 


D .   SUMMARY 

In  summary,  the  equations  used  in  the  analysis  are 


c2V2<{>  =    4>  (19) 


11] 


EI  — jr  w(x,t)    +   mw(x,t) 
9x 


-P 


(20) 


P   =   -    P<f> 


(21) 


w(x,t)    =   —   <t> 
8n 


(22) 
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III.   NUMERICAL  METHOD  OF  ANALYSIS  -  ACCELERATION  FORMULATION 

A.   FINITE  DIFFERENCING  SCHEME 
1.   The  Fluid 

The  wave  equation  is  hyperbolic  in  nature  and  for 
this  reason  an  explicit  formulation  for  the  timewise  response 
is  chosen.   The  conventional  central  finite  difference 
equations  are  used  in  both  space  and  time,  and  when  applied 
to  Eq.  (19)  give  [Ref.  3] 


2 
c 


(Ax) 


2   {*i,j-l  +  *i,j+l  "  4*i,j  +  *i+lj  +  ♦l-l,j} 


(AT)2   \*k-l  "  2<{,k  +  *k+l}  ij  (25) 

where  Ax  =  Ay  in  a  square  mesh,  i  and  j  are  mesh  point 
reference  indices,  and  k  is  the  time  step.   4>i  .  k+1  is  the 
only  unknown.   All  other  terms  are  known  or  can  be  evaluated 
prior  to  the  calculation  of  $.  .  . +,  . 
2.   The  Beam  Equation 

Conventional  central  finite  difference  schemes  are 
also  used  for  the  walls  in  both  space  and  time.   The  resulting 
equation  is  given  by  [Ref.  3] 
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-^  {w^  -  Uw1_1  +  6w±  -  4wi+1  +  w1+2J   + 


m 


(At) 


K-i  -  2wk +  wk+i}1  =  -pi, 


2  \"k-l     k    k+lj,      *Xk  (26) 


where 


Pi,k  "  "  2At   {~*k-l  +  ^k+l}.  . 


(27) 


and  j  =  1  or  j  =  Jmax  .   See  Fig.  1. 

A  similar  pair  of  equations  can  be  written  for  the  wall 
along  the  boundary  where  1=1.   In  Eq.  (26)   w.  ...  is  the 
only  unknown  in  the  explicit  scheme.   Note  that  the  pressure 
given  by  Eq .  (27)  is  a  function  of  the  fluid  potential  at 
the  wall. 

3-   The  Calculation  Scheme 

When  the  two  systems  are  combined,  the  other  linking 
mechanism  is  the  interface  condition  on  the  fluid  and  wall 
velocities.   Using  central  finite  differences,  Eq .  (22)  is 
given  by  [Ref.  3] 

2AT  fWi,k-l  +  Wi,k+l}  =  2AT  {"♦i.j-l.k  +  *ij+l,k}   (28) 
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In  the  calculation  scheme  the  wall  deflection  at  each  new 
step  k+1  is  computed  first,  and  then  the  velocity  potential 
is  computed  throughout  the  fluid.   In  order  to  compute  the 
wall  displacement  at  k+1  the  pressure  must  be  known  at  the 
preceding  time  step  k,  which  requires  knowledge  of  the 
velocity  potential  at  the  wall  at  the  new  time  step,  k+1, 
according  to  Eq.  (27).   The  solution  to  this  problem  is  to 
solve  for  the  wall  displacement  by  implicitly  solving  for 
the  new  velocity  potential  at  the  wall.   Thus,  using 
Eq .  (25)  to  eliminate  <j>.  .  .,,  from  Eq .  (27),  and  using 
Eq.  (28)  to  eliminate  <}>.  .  ,  .  (when  j  =  l)  from  the  result 


gives 


pi,k  =   "   2At   {"♦ij.k-l   "    ♦i,J,fc-l   +   2<f>i,j,k  + 


(tf)2  [*i-i,j,k  +  ♦i+i,j,k  •  Hij,k  +  *i,j+i,k 


*i,J+l,k  +  zf  (wi,k-l)]} 


A  similar  equation  can  be  obtained  for  the  wall  where 

j  =  j    .   Substituting  p  given  by  Eq .  (29)  into  Eq .  (26), 

fflclX 

combining  terms,  and  solving  for  w.  .+_  gives 
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wi,k+l 


■&-   At  \2$.    ,  ,  -  2(f).  ,  ,  .  + 
2m     ^  Ti,j,k    yi,j,k-l 


(TT}   c*i-i,j,k  +  *i+i,i,k  -  ^i.l.k  +  2*i,2,k 


+  -rr-    W 
At 


i,k^ 


+  2Wi,k  "  Wi,k-1 


EI  ,  At   v2  /         i.        ,  c 

(Ax)     I  i-2,k     l-l, k     l,k    i+I, k 


Wi+2,k} 


,   /,    ,   1   At2     2v  ,     , 

7  (1  +  S  A3T  c  }         (30) 


Equation  (30)  Is  written  for  the  wall  where  j  =  1,  0  <_  x  <  aH 
For  x  >  aH,  p  =  0.   Once  the  wall  deflection  is  known  at  k+1, 
the  interface  conditions  on  the  fluid  are  known.   The  calcula- 
tion of  the  fluid  velocity  potential  is  now  quite  straight- 
forward using  Eq.  (25)  and  Eq .  (28)  at  an  interface. 

On  the  free  surface  of  the  fluid  the  velocity 
potential  may  be  excited  by  an  applied  pressure.   In  this 
case,  applying  central  finite  differences  to  Eq.  (24)  gives 


*k+i  ■  *k-i  -  2  f  (pb>k  (31) 


where  (Pb)^  i-s  the  pressure  exerted  on  the  boundary  at  time 


step  k. 
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B.  INITIAL  CONDITIONS 

1.  The  Fluid 

Since  the  governing  partial  differential  equation 
(Eq.  (19))  is  second  order  with  respect  to  time,  two  initial 
conditions  are  required;  1)  the  velocity  potential  at  time 
zero",  and  2)  the  first  time  derivative  of  the  velocity 
potential,  which  is  proportional  to  any  initial  pressure  in 
the  fluid. 

2.  The  Wall 

The  beam  equation  (Eq.  (20))  is  also  second  order 
with  respect  to  time,  and  two  initial  conditions  are  required; 
1)  the  initial  displacement  of  the  wall,  and  2)  the  initial 
velocity  of  the  wall. 

C.  THE  TIME  STEP  FOR  NUMERICAL  STABILITY 
1.   The  Fluid 

For  hyperbolic  systems  solved  using  Eq .  (25)  3  the 
non-dimensional  time  step  ratio 


&  (32) 


where  primes  denote  non-dimensional  quantities,  must  be  less 
than  1  for  numerical  stability  for  a  one-dimensional  system 
[Ref .  3]  and  must  be  less  than   \J  2/2   for  numerical 
stability  for  a  two-dimensional  system  [Ref.  4]. 
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Non-dimensional  time  t'  is  given  by 


f  =  4r  (33) 


where  L  is  a  reference  length.   The  non-dimensional  length 
x'  is  given  by 


x'  =  r  (34) 


Substituting  Eqs .  (33)  and  (3*0  into  Eq .  (32)  gives 


*  <  -H=T  (35) 

c  vl  2 


for  numerical  stability  of  the  two-dimensional  system. 
2.   The  Wall 

Non-dimensionalizing  the  time  and  length  coordinates 
in  the  beam  equation  gives  [Ref.  3] 


t»  =  4-   \  51  (36) 

L2   II  m 


x»  -  g  (37) 

The  beam  equation  is  parabolic  in  nature  [Ref.  31  and 
requires  a  non-dimensional  time  step  ratio 


1   <  i  (38) 


(Ax') 


■  >2  -  2 
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for  numerical  stability.   Substituting  Eqs.  (36)  and  (37) 
into  Eq.  (38)  gives 


1  >   At 


2  "  (Ax)2 


S¥  (39) 


m 


Therefore 


At  <  (A^l   |g  (4o) 


3.   The  Combined  System 

The  assumption  has  been  made  that  the  fluid  will 
require  the  smallest  time  step  for  stability.   Thus,  there 
is  a  restriction  on  the  value  of  the  parameters  used  in  the 
combined  analysis,  in  particular  the  physical  parameters  of 
the  beam.   This  restriction  is  given  by 


<  4£  \\£r  (41) 


cIT  "  2 

according  to  Eqs.  (35)  and  (40) 
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IV.   NUMERICAL  METHOD  OF  ANALYSIS  -  VELOCITY  FORMULATION 

A.  "VELOCITY  FORMULATION  APPROACH 

Kreig  and  Monteith  [Ref.  5]  have  suggested  a  modification 
to  the  timewise  finite  differencing  scheme  outlined  in 
Section  III.   This  modification  is  described  as  a  "velocity" 
formulation. 

Consider  their  example  of  the  one-dimensional  spring  mass 
system. 

Mw  +  Kw   =   0  (42) 

Using  central  finite  differences,  Eq .  (42)  becomes 


M 


(At) 


2  K-l  "  2wk  +  Wk+l}  +  Kwk  =   °  (ZJ3) 


or 


"»i  ■  [2  - 1  (At,2]wk  -  -Vi  <M) 


Now  consider  the  velocity  formulation 


Mv  +  Kvj   =   0 


W   =   V 


(15) 
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Defining  the  velocity  at  half  time  steps  and  the  displacement 
at  whole  time  steps  leads  to 


At 


{"  v.*,  +  VbJ  +  KWk  =  °  •  W) 


or 


vk+i,  =  Vk-H   "  I  At  wk  ^7) 


and 

wk+l  =  wk  +  At  vk+%  (J48) 

With  either  formulation  the  time  step  is  the  same,  and  in 
the  real  number  system  without  round-off  error  the  results 
would  be  identical.   However,  computers  do  not  compute  in 
the  real  number  system.   Round-off  errors  are  always  present 
Consider  the  term  in  brackets  in  Eq .  (44) 


[2  -  |  (At)2] 


Depending  on  the  number  of  digits  retained,  the  significance 

K      P 
of  jjr  (At)   could  be  completely  lost  in  round  off. 


24 


B.   FINITE  DIFFERENCING  SCHEME 
1.   The  Fluid 

The  fluid  equation  derived  in  Section  III  can  be 
written  as 


c2V2<J>  =  41  (i»9) 


♦  ■--•-*  (50) 

A  central  finite  difference  scheme  is  proposed  centering 

the  "velocity"  at  whole  time  steps  and  the  "displacement" 
at  half  time  steps.   Note,  the  pressure  is  directly  computed 
at  each  whole  time  step  since 

Pk  =  -  P  ^k  (51) 

The  finite  difference  equation  for  the  fluid  is  given  by 


2 

c 


(Ax 


^2   {♦l,J-l  +  *i,J+l  "  "♦l.J  +  *i+l9j  +  ♦i-l,J/kji  = 
AT   K-l  +  \}    .  ,  (52) 


where  \\> .    .  ,  is  the  only  unknown.   The  potential  is  obtained 
from  the  finite  difference  expression  for  Eq .  (50) 


♦ij.k+ii  =^i,J,k-%  +  At*ij,k  (53) 
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and 


2.   The  Wall 

Similarly,  the  beam  equation  can  be  written 

EI  — jj-  w(x,t)  +  mv  =  -p  (54) 

dx 

w(x,t)  =  v  (55) 


For  these  equations  a  central  finite  difference  scheme  is 
proposed  centering  the  velocity  v  at  half  time  steps  and 
the  displacement  w  at  whole  time  steps .   Note  that  the 
velocity  needed  for  the  interface  condition  is  directly 
computed, 

Vi,k^  =  U   {-♦l.J-1  +  *l,J+l}k^  (56) 

The  finite  difference  equation  for  the  beam  is 


(57) 


it  <w.  „  -  4w.  ,  +  6w.  -  4wj.n  +  w.  ,„> 

,.     vl    1-2      l-l     i      l+l     i+2], 
(Ax)   l  Jk 


AT  {-vi,k-%  +  Vi,k+J   =   p  *ifj,] 


The  unknown  is  v.  ,  ,,  .   From  Eq .  (55) 

1  ,K  +  *2 


wi,k+l  =  wi,k  +  At  vi,k+^  (58) 
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3.   The  Combined  System 

With  the  velocity  function  there  is  no  need  to  solve 
for  the  beam  displacement  in  terms  of  an  Implicit  solution 
of  the  potential  at  the  wall.   In  the  velocity  formulation 
the  pressure  is  computed  directly;  no  differencing  is 
required.   Figure  2  graphically  illustrates  the  time  step 
centering  of  both  the  acceleration  formulation  and  the 
velocity  formulation. 

The  solution  scheme  for  the  velocity  formulation  is: 

1.  Solve  for  the  wall  velocity,  Eq .  (57),  based  on  the 
known  pressure  at  the  wall  and  the  previous  displace- 
ment and  velocity. 

2.  Calculate  the  new  wall  displacement,  Eq .  (58),  based 
on  the  new  wall  velocity  and  the  previous  displacement. 

3.  Solve  for  ty}   Eq.  (52),  based  on  the  previous  tya    the 
new  wall  velocity  (for  the  interface  conditions)  and 
the  previous  (j>. 

4.  Calculate  the  new  <J>,  Eq .  (53),  based  on  the  previous 
$  and  the  newly  calculated  lfi. 

5-   Return  to  step  1  and  repeat. 
As  in  the  example  given  by  Kreig  and  Monteith  this  formula- 
tion is  exactly  equal  to  the  acceleration  formulation  if  the 
beam  and  fluid  equations  are  considered  separately.   However, 
in  the  combined  system  the  fluid  equation  is  displaced 
exactly  one-half  time  step  from  the  fluid  equation  used  in 
the  acceleration  formulation.   In  the  acceleration  formulation 


27 


ACCELERATION  FORMULATION 


4> 


-1 


w 


0 


+1 


+2 


*    TIME 


+3 


£ 


W 


VELOCITY   FORMULATION 


-1 


0  C$\  $ 


\  <p   <f>\  $($ 


\ 


+1 


1 


+2 


+3 


T 


-V    TIME 


FIG.    2 


23 


the  pressure  depends  on  (J>k+1  and  <f>k_-,  v;hereas  in  the  velocity 
formulation,  it  depends  on  <{>    and  <(>,    .   Thus,  even  with 
no  round-off  error  there  will  be  a  difference  between  the  two 
solutions.   Refer  to  Part  D  for  a  discussion  of  the  signif- 
icance of  this  difference. 

C.   INITIAL  CONDITIONS 

As  with  the  acceleration  formulation  a  total  of  four 

initial  conditions  are  required.   However,  recall  that  two 

of  the  four  variables  are  defined  one-half  time  step  on 

either  side  of  t  =  0.   They  are  the  wall  velocity  and  the 

velocity  potential  of  $   in  the  fluid. 

The  actual  initial  conditions  can  be  denoted  as  w  and 

o 

<j)   .   The  computer  code  assumes  that 

v  .  =  wrt  (59) 


-%  "O 


and 


«    =  0o  (60) 


D.   THE  TIME  STEP  FOR  NUMERICAL  STABILITY 

Theoretically,  since  the  same  parameters  are  used  in 
both  formulations,  the  time  step  for  numerical  stability 
should  be  exactly  the  same  for  both  formulations.   In  fact, 
this  was  found  to  be  true  for  the  beam  equation  and  for  the 
fluid  equation  when  each  was  tested  separately.   However, 
when  the  systems  were  coupled,  and. the  maximum  theoretical 
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time  step  was  used,  substantial  instability  exhibited 

itself  immediately  in  the  computer  results.   The  instability 

occurred  at  the  interface  and  it  occurred  within  two  time 

steps  of  a  pressure  pulse  striking  a  wall.   A  subsequent 

velocity  was  imparted  to  the  wall  which  was  always  large 

enough  to  cause  a  reversal  in  the  fluid  pressure  at  the  wall 

at  the  next  time  step. 

A  primitive  but  effective  solution  to  the  instability 

problem  was  obtained  by  reducing  the  time  step  arbitrarily 

to 

At  =   A*  (61) 

c  i   3 
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V.   DESCRIPTION  OF  THE  COMPUTER  PROGRAMS 

Two  Fortran  IV  programs  were  written  and  tested.   The 
theory  behind  each  is  outlined  in  Sections  III  and  IV  above. 
The  two  programs  are  quite  similar  in  construction  and 
require  similar  amounts  of  core  storage.   C.P.U.  time 
required  is  also  similar.   Only  one  program  will  be  described 
in  detail,  but  differences  between  the  two  will  be  pointed 
out  where  they  occur. 

A.   GENERAL  FEATURES 

The  programs  allow  for  complete  flexibility  in  any  of 
the  possible  initial  conditions,  both  on  the  walls  and  on 
the  fluid.   At  present  the  program  allows  for  three  walls 
and  a  variable  fluid  depth.   The  ends  of  each  wall  are 
considered  to  be  clamped.   The  wall  boundary  conditions  can 
be  changed  by  any  future  user  by  changing  six  equations  in 
the  velocity  formulation  and  twelve  in  the  acceleration 
formulation.   The  location  of  these  equations  will  be  pointed 
out  below.   The  program  provides  for  a  non-dimensional  thirty 
by  thirty  mesh  point  domain,  with  the  reference  length  being 
the  length  of  the  horizontal  wall  (See  Fig.  1).   All  variables, 
vectors  and  arrays  are  placed  in  blank  Common,  to  allow  for 
complete  flexibility  in  communication  between  subroutines. 
Further  all  vectors,  and  arrays  are  initialized  to  zero  value 
before  analysis  begins. 
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Subroutines  are  provided  to: 

1.  Compute  constants  derived  from  input  parameters, 
normalize  the  dimensions  of  the  problem,  and  compute 
or  set  initial  conditions. 

2.  Output  data  obtained  for  each  mesh  point  in  the  domain. 
The  format  is  tabular. 

3.  Compute  the  beam  equations. 

*1 .   Compute  the  fluid  equations. 
5.   Provide  for  a  special  graphical  output. 
Each  subroutine  will  be  described  in  detail. 


B.   MAIN  PROGRAM 

1.   Variable  Description  —  Velocity  Formulation 


wa(i; 

WB(I; 

wc(i: 
va(i: 

VB(I 

vc(i; 
ua(i; 
ub(i; 

UC(I 

XA(I 

XB(I 

XC(I 

WADOT(I) 

WBDOT ( I ) 


-  left  vertical  wall  displacement  at  k 

-  horizontal  wall  displacement  at  k 

-  right  vertical  wall  displacement  at  k 

-  left  vertical  wall  displacement  at  k+1 

-  horizontal  wall  displacement  at  k+1 

-  right  vertical  wall  displacement  at  k+1 

-  left  vertical  wall  velocity  at  k-Jj 

-  horizontal  wall  velocity  at  k-h 

-  right  vertical  wall  velocity  at  k-% 

-  left  vertical  wall  velocity  at  k+% 

-  horizontal  wall  velocity  at  k+h 

-  right  vertical  wall  velocity  at  k+% 

-  initial  left  vertical  wall  velocity 

-  initial  horizontal  wall  velocity 
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WCDOT(I) 

PRA(I) 

PRB(I) 

PRC(I) 

PK(I) 

IWA(I) 

IWB(I) 

IWC(I) 


-  initial  right  vertical  wall  velocity 

-  initial  left  vertical  wall  pressure 

-  initial  horizontal  wall  pressure 

-  initial  right  vertical  wall  pressure 

-  fluid  free  surface  applied  pressure 

-  integer  value  of  left  vertical  wall 
displacement 

-  integer  value  of  horizontal  wall 
displacement 

-  integer  value  of  right  vertical  wall 
displacement 


P(I,J) 

PN(I,J) 
PL(I,J) 

PRESS(I,J) 

PHIDOT(I,J) 

IPRESS(I.J) 


-  value  of  <J>  in  the  fluid  at  k-% 

-  value  of  <j>  in  the  fluid  at  k+% 

-  value  of  ty   in  the  fluid  at  k-1 

-  value  of  pressure  in  fluid  at  k 

-  initial  value  of  \p   in  the  fluid 

-  integer  value  of  pressure  in  the  fluid. 

In  the  displacement  formulation  the  following  variables  differ, 
VA,VB}VC        are  not  used 


WA,WB,WC 

XA,XB,XC 

UA,UB,UC 

P(I,J) 

PN(I,J) 
PL(I,J) 


are  values  of  wall  displacement  at  k 
are  values  of  wall  displacement  at  k+1 
are  values  of  wall  displacement  at  k-1 

are  values  of  $   at  k 
are  values  of  $  at  k+1 
are  values  of  <J>  at  k-1 
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2.   Main  Program  Functions 

In  the  velocity  formulation  the  main  program 
accomplishes  the  following  functions: 

1.  Controls  calls  to  other  subroutines. 

2.  Computes  first  value  of  PN(I,J)  from  initial  conditions. 

3.  Saves  the  value  of  the  displacement  at  a  selected 
wall  mesh  point  at  selected  time  steps  for  later 
plotting  as  a  time  vs.  displacement  graph. 

4.  Turns  off  a  free  surface  pressure  disturbance  at  any 
preselected  time  step. 

In  addition  the  "acceleration"  formulation  program 

5.  Computes  the  k  =  1  value  of  wall  displacement  from 
initial  conditions. 

C.   SUBROUTINE  INP(IER) 
This  subroutine 

1.  Reads  all  input  data. 

2.  Computes  constants. 

3.  Assigns  values  to  initial  conditions  other  than  zero 
value . 

1 .   Variable  Description 

RHOF  -  Density  of  the  fluid 

E  -  Young's  Modulus 

RL  -  Horizontal  length 

RH  -  Vertical  height  where  RH  <_  RL 

AL  -  Fractional  depth  of  fluid  0  <  AL  <  1.0 

C  -  Speed  of  sound  in  the  fluid. 
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TH 
RHO 

NH1 
Nl 

DELX 
DELT 
RIO 

NA 

C1,C2,C3 
04,05,06,07,08,09,0 

D1,D2,D3 

NH 

NHI 

N 

NAI 


-  Thickness  of  beam 

-  Density  of  beam  material  in  lbs 
per  in3 

-  Number  of  vertical  mesh  points 

-  Number  of  horizontal  mesh  points 
where  NHI  <  Nl  =  30 

-  Width  of  mesh 

-  Time  step 

-  Computed  value  of  I,  moment  of 
inertia 

-  Number  of  vertical  mesh  points  in 
the  fluid. 

-  Constants  common  to  both  programs 

-  Constants  used  only  in  acceleration 
formulation 

-  Constants  used  in  velocity  formulation 

-  Equals  NHI  -  1 

-  Equals  NHI  -  2 

-  Equals  Nl  -  1 

-  Equals  NA  -  1 


2 .  Input  Cards 

Only  two  input  cards  are  used  to  input  RL,  RH,  AL, 
C,  TH,  RHO,  RH0P,  E.   The  format  will  be  described  in 
Appendix  A. 

3.  Non-Zero  Initial  Conditions 

Any  non-zero  initial  conditions  are  inserted  by  the 
user  at  the  end  of  Subroutine  INP . 
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D.  SUBROUTINE  OUTP(IER) 

Subroutine  OUTP  prints  out  in  tabular  format  the  values 
of  each  wall  mesh  point  displacement  and  the  pressure  at 
each  mesh  point  throughout  the  domain. 

E.  SUBROUTINE  WALL(IER) 

Subroutine  WALL  computes  the  values  of  the  wall  velocity 
and  displacement  in  the  velocity  formulation  and  the  wall 
displacement  in  the  acceleration  formulation.   If  a  future 
user  desired  to  change  the  boundary  conditions  of  the  beam 
equation,  then  any  equation  with  a  computed  answer  stored 
in  a  variable  indexed  with  a  2,    N,  or  NH  would  have  to  be 
changed.   In  the  velocity  formulation  only  the  equations 
for  VA(2),  VA(NH),  VB(2)5  VB(N) ,  VC(2),  VC(NH)  need  to 
be  changed.   In  the  acceleration  formulation  these  equations 
are  for  XA(2),  XA(NH),  XB(2),  XB(N),  XC(2),  XC(NH)  and  recall 
that  these  equations  must  also  be  changed  in  the  main  program, 

F.  SUBROUTINE  FLUID(IER) 

Subroutine  FLUID  computes  the  value  of  ty   and  <f>  at  k  and 
k+h   respectively  in  the  velocity  formulation  and  computes 
the  values  of  <J>  at  k+1  in  the  acceleration  formulation. 
Pressure  at  each  mesh  point  is  also  computed  in  FLUID. 

G.  SUBROUTINE  GRAPHO(IER) 

Subroutine  GRAPHO  is  an  attempt  to  make  some  sense  out 
of  the  possible  1900  numerical  values  calculated  at  each 
time  step.   The  subroutine  locates  an  absolute  maximum  value 
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for  wall  displacement  and  fluid  pressure  and  from  these 
values  normalizes  the  entire  domain  to  integer  values  from 
-99  to  +100.   This  information  is  then  presented  in  a  format 
quite  similar  to  Fig.  (1) .   An  integer  value  is  placed  at 
each  mesh  point  in  the  domain. 
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VI.   SOLUTIONS 

A.   CHECK  CASES 

In  order  to  determine  if  the  programs  correctly  modeled 
the  system  described  above,  several  test  cases  were  run  on 
each  formulation. 
1.   The  Fluid 

By  controlling  the  subroutines  called  by  the  main 
program,  the  user  is  able  selectively  to  compute  results 
for  the  fluid  alone  or  the  wall  alone.   One  test  case  for 
the  fluid  alone  was  a  step  pressure  disturbance  applied 
uniformly  over  the  entire  upper  surface  for  a  finite  time 
interval.   The  walls  were  treated  as  rigid.   This  test 
showed  that  the  pulse  propagation  and  reflection  behaved 
as  expected  in  both  formulations.   The  uniform  pulse  traveled 
uniformly  through  the  fluid.   The  value  of  pressure  doubled 
as  expected  at  the  horizontal  rigid  wall  and  returned 
uniformly  to  the  free  surface.   See  Figs.  3  and  k. 

Another  test  case  consisting  of  a  point  pressure 
disturbance  at  the  middle  of  the  free  surface  showed  that 
two-dimensional  pulse  propagation  in  the  fluid  was  as 
expected.   The  amplitude  of  the  pressure  decreased  radially 
from  the  point  source. 

In  both  test  cases  the  theoretical  time  step  for 
stability  proved  to  be  valid  in  both  formulations. 
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2.   The  Wall 

A  free  vibration  analysis  was  conducted  using  the 
beam  equation  to  test  the  wall  response  alone.   The  wall 
was  initially  displaced  to  conform  to  the  first  mode  of 
vibration.   Rogers  [Ref.  1]  gives  a  value  of  51.3  cycles  per 
second  for  the  value  of  the  frequency  of  the  first  mode  of 
vibration  for  the  beam  parameters  used  in  the  test.   The 
results  from  both  formulations  returned  a  first  mode  of  vibra- 
tion frequency  of  approximately  50  cycles  per  second.   Again 
the  theoretical  time  step  for  stability  proved  to  be  valid. 

B.   EXAMPLE  PROBLEMS 

The  types  of  problems  these  two  codes  were  developed  to 
solve  involve  the  combined  system.   The  first  example  prob- 
lem consisted  of  a  uniform  pressure  pulse  at  the  free  sur- 
face of  a  flexible  tank  for  2.26  milleseconds .   When  a  solu- 
tion was  first  attempted  using  the  velocity  formulation,  the 
previously  mentioned  instability  was  discovered  when  the 
theoretical  time  step  for  numerical  stability  of  the  fluid 
was  used.   Since  no  other  solutions  to  this  type  of  problem 
were  available,  first  priority  centered  on  determining  the 
cause  and  cure  of  the  instability.   A  reduction  in  the  time 
step  as  discussed  in  Section  IV  solved  the  instability  prob- 
lem.  The  velocity  formulation  was  run  with  a  reduced  time 
step,  and  the  acceleration  formulation  was  run  using  the 
theoretical  time  step.   The  results  for  the  displacement 
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of  the  center  of  the  horizontal  wall  are  shown  in  Fig.  (5). 
The  dominant  frequency  of  vibration  is  approximately  12 
cycles  per  second  for  both  formulations.   There  is  a  small 
amplitude  discrepency  of  about  5%   between  the  two  formula- 
tions.  The  resulting  vibration  appeared  to  involve  all 
modes  in  the  beam. 

The  acceleration  formulation  was  also  used  to  solve  the 
same  problem  with  a  uniform  pulse  held  for  .226  milleseconds 
or  1/10   the  previous  example.   Figure  (5)  again  denotes 
a  frequency  of  approximately  12  cycles  per  second,  and  the 
amplitude  is  proportional  to  the  previous  case.   Maximum 
amplitude  is  about  one-tenth  that  of  the  previous  case. 

It  is  interesting  to  note  that  if  the  fluid  is  considered 
as  incompressible,  and  its  mass  is  added  to  the  mass  of  the 
beam,  the  frequency  of  free  vibration  as  given  by  Rogers 
[Ref.  2]  is  11.1  cycles  per  second. 

Another  example  problem  was  run  with  a  pressure  exerted 
at  the  two  center  mesh  points  of  the  free  surface  for  approx- 
imately 450  microseconds.   In  this  case,  the  first  mode  of 
vibration  predominated  in  the  horizontal  wall  at  a  frequency 
of  13  cycles  per  second.   Again  only  minimal  differences 
in  amplitude  were  observed  between  the  two  formulations. 
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VII.   CONCLUSIONS 

The  results  obtained  to  date  appear  to  support  the 
following  views: 

1.  The  amplitude  of  vibration  in  the  walls  are  depen- 
dent upon  the  length  of  time  the  free  surface  excitation 
pressure  is  applied.   This  would  be  expected  purely  from 
energy  concepts . 

2.  Choice  of  formulation  is  up  to  the  user.   Each  has 
its  own  virtues  and  its  vices.   The  velocity  formulation 
requires  less  core  storage  and  less  CPU  time  for  the  same 
number  of  time  steps.   (The  acceleration  formulation  requires 
15%   more  CPU  time  than  the  velocity  formulation;  however, 

the  time  step  required  is  19  percent  smaller  in  the  velocity 
formuation.)   The  velocity  formulation  would  be  easier  to 
modify . 

3.  No  claims  are  made  as  to  the  efficiency  of  either 
program  as  presently  written.   In  fact  it  is  estimated  that 
the  program  size  could  be  reduced  considerably. 

ty .   More  work  needs  to  be  done  in  checking  out  other 
simple  cases  to  prove  the  validity  of  the  codes,  for  exam- 
ple, pressure  pulses  passing  through  the  fluid  at  angles 
from  the  horizontal. 

5.   Either  program  will  give  a  future  user  valuable 
insight  about  structural  behavior  in  contact  with  a  fluid. 


Hl\ 


APPENDIX  A 
USER'S  INSTRUCTIONS 

Input  of  Initial  conditions  has  already  been  discussed 
in  Section  V  above.  Actual  input  data  is  contained  in  two 
cards  with  the  following  format : 


Card  1 

Column 

Variable 

Card  1: 

1-10 

RL 

11-20 

RH 

21-30 

AL 

Definition  Format 

Reference  length       F10.4 

Vertical  height        F10.4 

Fractional  depth       F10.4 
of  fluid 

31-40  C  Velocity  of  sound      F10.4 

in  the  fluid 

41-50  TH         Thickness  of  wall      F10.4 

51-60  RHO        Density  of  wall  _      F10.4 

material  in  #/in 

Card  2: 


1-12  RHOF       Density  of  fluid       E12.4 

in  #-sec2/in^ 

13-24  E  Young's  Modulus  of     E12.4 

wall  material 


45 
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